This is the first of two Required Practice activities for today and this one that focuses on introducing the plotly package. In my opinion, the plotly package is most useful for generating charts (e.g., scatterplots) for inclusion in a dashboard or webpage, but you can also use it generate maps too. It is just waaaaaay easier to produce interactive maps with the tmap package 😉


There is often value in working with a familiar dataset in times like this. You have a shapefile of (redacted) flood insurance claims paid in the state of Virginia compiled by the Federal Emergency Management Administration included in the .zip file. You can read more about the source of these data here if you like but we saw this dataset last week as part of the exam review.


Below, we read in the shapefile from hard disk using the sf::st_read function…


va_nfip_claims_0 <- st_read("./NFIP_Claims_VA.shp")
## Reading layer `NFIP_Claims_VA' from data source 
##   `C:\Users\bw6xs\Documents\PLAN 6122\2023\Lectures\Week 9\Week 9 - Lecture 1 - Required Practice\NFIP_Claims_VA.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2959 features and 38 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -82.8 ymin: 36.5 xmax: -75.3 ymax: 39.4
## Geodetic CRS:  WGS 84


Now let’s practice generating some common graphics using plotly. The primary advantage here is that these graphics are interactive by default because this R package allows us to leverage functionality from the plotly JavaScript library commonly used for web development.


va_nfip_claims_0  %>%
  plot_ly() %>% 
  add_histogram(x = ~yrOfLss)
va_nfip_claims_0  %>%
  filter(yrOfLss >= 2003) %>%
  plot_ly() %>% 
  add_histogram(x = ~yrOfLss, marker = list(color = "orange", line = list(color = "royalblue", width = 2))) %>%
  layout(title = "Flood Insurance Claims in VA", plot_bgcolor = "grey20", 
         xaxis = list(title = "Year"),
         yaxis = list(title = "Number"))


Note that we can chain plotly code together with dplyr functions for data wrangling in the same way that we saw with ggplot2 graphics 💪. In the above chunk, the plot_ly function initializes a plotly object and the add_histogram function builds a histogram based on the number of flood insurance claims filed in Virginia by year. The second example just shows how to tailor the resulting histogram chart. Note that the marker and line arguments take a list with multiple parameters and that both resulting graphics are dynamic in that we get additional information when we mouse over (hover) the bars.


# Boxplot...
va_nfip_claims_0 %>%
  filter(occpncT < 4) %>%
plot_ly() %>%
  add_boxplot(x = ~occpncT, y = ~amnPOBC, 
              text = ~paste(occpncT, floodZn),
              hoverinfo = "text") %>%
  layout(title = "Flood Insurance Claims Paid for Building Costs", plot_bgcolor = "grey20", 
         xaxis = list(title = " ", 
                      ticktext = list("Single Family Res", "2 to 4 Unit Res", "More Than 4 Unit Res"), 
                      tickvals = list(1, 2, 3)), 
         yaxis = list(title = "Amount Paid ($)"))
# Scatterplot...
va_nfip_claims_0  %>%
  drop_na(amnPOBC) %>%
  filter(ttlBlIC <= 250000) %>%
  plot_ly() %>%
  add_trace(y = ~amnPOBC, x = ~ttlBlIC, marker = list(color = "red"),
            text = ~paste("My policy was for", scales::dollar(ttlBlIC), "but I received ", scales::dollar(amnPOBC), "after the storm." ),
            hoverinfo = "text") %>%
  layout(title = "Flood Insurance Claims Paid for Building Costs", plot_bgcolor = "grey20", 
         xaxis = list(title = "Total Insurance Amount"), 
         yaxis = list(title = "Amount Paid Bldg Claim"))


If we want to access cool basemaps—and of course we do—it is best to sign up for MapBox token by visiting this page.




Sys.setenv('MAPBOX_TOKEN' = 'PASTE_YOUR_STRING_HERE')

Copy the line of code shown above into the Console, insert your token from the MapBox website, then press Enter to save it into memory BEFORE running the code chunk below.


What do you call a mouse that swears?

A cursor. 😆

# Dynamic map with basemap included...
va_nfip_claims_0 %>%
  drop_na(ttlBlIC, amnPOBC) %>%
plot_mapbox(color = ~amnPOBC, 
            text = ~paste("My policy was for", scales::dollar(ttlBlIC), "but I received ", scales::dollar(amnPOBC), "after the storm." ),
            hoverinfo = "text") %>%
  layout(mapbox = list(style = 'carto-positron', zoom = 10, center = list(lat = ~median(st_coordinates(va_nfip_claims_0)[, 2]), 
                                                                             lon = ~median(st_coordinates(va_nfip_claims_0)[, 1]))),
            title = " ") %>% 
  colorbar(title = "Amount Paid ($) <br> Building Claim")


Now let’s demonstrate how to coerce a ggplot map to a plotly map… 😉


# Apply a map projection... maybe UTM Zone 17 North...
va_nfip_claims_sf <- st_transform(va_nfip_claims_0, crs = "EPSG:32617")

# Filter out records where year of loss is 2003
va_nfip_claims_2003_sf <- va_nfip_claims_sf %>%
  filter(yrOfLss == 2003)


# Get a layer with boundaries for context...
bound_0 <- st_read("https://opendata.arcgis.com/datasets/e3c8822a4adc4fc1a542a233893a46d4_0.geojson")
## Reading layer `SDE_USDC_CENSUS_VA_COUNTY' from data source 
##   `https://opendata.arcgis.com/datasets/e3c8822a4adc4fc1a542a233893a46d4_0.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 134 features and 20 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -83.67539 ymin: 36.54076 xmax: -75.16644 ymax: 39.46601
## Geodetic CRS:  WGS 84
# Apply a map projection... maybe UTM Zone 17 North... to be consistent!
bound <- st_transform(bound_0, crs = "EPSG:32617")

# Map claims by structure type AND payment amount, KNOWING that the vast
# majority of structure were single-family residential...
claims_2003 <- ggplot() + 
  geom_sf(data = va_nfip_claims_2003_sf, aes(color = as_factor(occpncT), size = amnPOBC), shape = 19, alpha = 0.4) +
  geom_sf(data = bound, color = "dodgerblue", fill = NA) +
  scale_color_viridis(discrete = TRUE, name = "Structure Type", labels = c("Single Family Res", 
                                                                           "2 to 4 Unit Res", 
                                                                           "More Than 4 Unit Res", 
                                                                           "Non-residential")) +
  scale_size_continuous(name = "Building Claim Paid($)") +
  coord_sf(crs = st_crs(4326)) +
  theme_void() + 
  theme(legend.position="bottom") 

ggplotly(claims_2003)


We don’t have that much control over the details of the plot if we go this route. Let’s create a map from scratch instead… 😝


# Group and count by county (locality)
by_county <- va_nfip_claims_2003_sf %>%
  drop_na(amnPOBC) %>%
  group_by(contyCd) %>%
  count()


# Join to existing sf object for plotting...
by_county_joined <- st_join(bound, by_county, by = c("GEOID" = "contyCd"), left = TRUE)

# Generate the map...
plot_ly(data = bound, stroke = I("grey"), showlegend = FALSE) %>%
    add_sf(data = by_county_joined, 
         split = ~NAME,
         color = ~n,
         colors = "Reds",
         stroke = I("grey"),
         text = ~paste0(NAME, "\n", n),
         hoveron = "fills",
         hoverinfo = "text", 
         inherit = FALSE) %>%
    layout(showlegend = FALSE) %>%
  colorbar(title = "Number of \n Claims")

Your Turn


Insert a new code chunk that generates one or two graphics of your choice using plotly. This is can be with the dataset we have here OR another dataset of your choosing.



You have reached the end! 👅